using System;
using MathLib;

namespace UraniaLib
{
	/// <summary>
	/// Summary description for UraniaPlanet.
	/// </summary>
	public class UraniaPlanet
	{
	
		public UraniaPlanet()
		{
			//
			// TODO: Add constructor logic here
			//
		}

		public static void CalcObjectPosFromSun(DateTime dEpoch, DateTime dDate, double fOrbPeriod, double fEcc, double fEpochLong, double fPerihelionLong, double fSMA, ref double fLong, ref double fRadius)
		{
			double fM;
			double fV;
			double fD;
			fD = UraniaTime.GetDaysBetween(dDate, dEpoch);

			fM = ((360.0/365.242191) * (fD/fOrbPeriod)) + fEpochLong - fPerihelionLong;
			fV = fM + ((360.0/Math.PI) * fEcc * Trig.Sin(fM));

			fLong = fV + fPerihelionLong;
			fRadius = (fSMA * (1 - (fEcc * fEcc)))/(1 + (fEcc * Trig.Cos(fV)));

		}

		public static void CalcPlanetCoord(DateTime dEpoch, DateTime dDate, double fPOrbPeriod, double fPEcc, double fPEpochLong, double fPPeriLong, double fPSMA, double fPAscNode, double fPIncl, double fEOrbPeriod, double fEEcc, double fEEpochLong, double fEPeriLong, double fESMA, bool bInnerPlanet,  ref double fRA, ref double fDecl)
		{
			double fD, fPN, fPM, fPL, fPV, fPR;
			double fEN, fEM, fEL, fEV, fER;
			double fPsi, fX, fY, fA, fL1, fR1, fLambda, fBeta, fT;

			fD = UraniaTime.GetDaysBetween(dDate, dEpoch);
			fPN = (360.0/365.242191) * (fD / fPOrbPeriod);
			fPN = Trig.PutIn360Deg(fPN);
			fPM = fPN + fPEpochLong - fPPeriLong;
			fPL = fPN + ((360.0 / Math.PI) * fPEcc * Math.Sin(Trig.DegToRad(fPM))) + fPEpochLong;
			fPL = Trig.PutIn360Deg(fPL);
			fPV = fPL - fPPeriLong;
			fPR = (fPSMA * (1 - (fPEcc * fPEcc))) / (1 + (fPEcc * Math.Cos(Trig.DegToRad(fPV))));

			fEN = (360.0/365.242191) * (fD / fEOrbPeriod);
			fEN = Trig.PutIn360Deg(fEN);
			fEM = fEN + fEEpochLong - fEPeriLong;
			fEL = fEN + ((360.0 / Math.PI) * fEEcc * Math.Sin(Trig.DegToRad(fEM))) + fEEpochLong;
			fEL = Trig.PutIn360Deg(fEL);
			fEV = fEL - fEPeriLong;
			fER = ((1 - (fEEcc * fEEcc))) / (1 + (fEEcc * Math.Cos(Trig.DegToRad(fEV))));

			fPsi = Math.Asin(Math.Sin(Trig.DegToRad(fPL - fPAscNode)) * Math.Sin(Trig.DegToRad(fPIncl)));
			fPsi = Trig.RadToDeg(fPsi);

			fY = Math.Sin(Trig.DegToRad(fPL - fPAscNode));
			fX = Math.Cos(Trig.DegToRad(fPL - fPAscNode));
	//		fT = Math.Atan((fY/fX) * Math.Sin(Trig.DegToRad(fPIncl)));
			fT = Math.Atan((fY/fX));
			fT = Trig.RadToDeg(fT) ;
			fT = Trig.TanQuadrant(fX, fY, fT);
			fL1 = fT + fPAscNode;
			fL1 = Trig.PutIn360Deg(fL1);
			fR1 = fPR * Math.Cos(Trig.DegToRad(fPsi));

			if (bInnerPlanet == true)
			{
				fA = Math.Atan((fR1 * Math.Sin(Trig.DegToRad(fEL - fL1)))/(fER - (fR1 * Math.Cos(Trig.DegToRad(fEL - fL1)))));
				fA = Trig.RadToDeg(fA);
				fLambda = 180 + fEL + fA;
				fBeta = Math.Atan((fR1 * Math.Tan(Trig.DegToRad(fPsi))* Math.Sin(Trig.DegToRad(fLambda - fL1)))/(fER * Math.Sin(Trig.DegToRad(fL1 - fEL))));
				fBeta = Trig.RadToDeg(fBeta);
			}
			else
			{
				fLambda = Math.Atan((fER * Math.Sin(Trig.DegToRad(fL1 - fEL)))/(fR1 - (fER * Math.Cos(Trig.DegToRad(fL1 - fEL)))));
				fLambda = Trig.RadToDeg(fLambda) + fL1;
				fLambda = Trig.PutIn360Deg(fLambda);
				
				fBeta = Math.Atan((fR1 * Math.Tan(Trig.DegToRad(fPsi))* Math.Sin(Trig.DegToRad(fLambda - fL1)))/(fER * Math.Sin(Trig.DegToRad(fL1 - fEL))));
				fBeta = Trig.RadToDeg(fBeta);
			}

			UraniaCoord.ConvEclToEqu(23.441884,fLambda, fBeta, ref fRA, ref fDecl);
		}

		public static void CalcPlanetPhase(DateTime dEpoch, DateTime dDate, double fPOrbPeriod, double fPEcc, double fPEpochLong, double fPPeriLong, double fPSMA, double fPAscNode, double fPIncl, double fEOrbPeriod, double fEEcc, double fEEpochLong, double fEPeriLong, double fESMA, bool bInnerPlanet,  ref double fPhase)
		{
			double fD, fPN, fPM, fPL, fPV, fPR;
			double fEN, fEM, fEL, fEV, fER;
			double fPsi, fX, fY, fA, fL1, fR1, fLambda, fT;

			fD = UraniaTime.GetDaysBetween(dDate, dEpoch);
			fPN = (360.0/365.242191) * (fD / fPOrbPeriod);
			fPN = Trig.PutIn360Deg(fPN);
			fPM = fPN + fPEpochLong - fPPeriLong;
			fPL = fPN + ((360.0 / Math.PI) * fPEcc * Math.Sin(Trig.DegToRad(fPM))) + fPEpochLong;
			fPL = Trig.PutIn360Deg(fPL);
			fPV = fPL - fPPeriLong;
			fPR = (fPSMA * (1 - (fPEcc * fPEcc))) / (1 + (fPEcc * Math.Cos(Trig.DegToRad(fPV))));

			fEN = (360.0/365.242191) * (fD / fEOrbPeriod);
			fEN = Trig.PutIn360Deg(fEN);
			fEM = fEN + fEEpochLong - fEPeriLong;
			fEL = fEN + ((360.0 / Math.PI) * fEEcc * Math.Sin(Trig.DegToRad(fEM))) + fEEpochLong;
			fEL = Trig.PutIn360Deg(fEL);
			fEV = fEL - fEPeriLong;
			fER = ((1 - (fEEcc * fEEcc))) / (1 + (fEEcc * Math.Cos(Trig.DegToRad(fEV))));

			fPsi = Math.Asin(Math.Sin(Trig.DegToRad(fPL - fPAscNode)) * Math.Sin(Trig.DegToRad(fPIncl)));
			fPsi = Trig.RadToDeg(fPsi);

			fY = Math.Sin(Trig.DegToRad(fPL - fPAscNode));
			fX = Math.Cos(Trig.DegToRad(fPL - fPAscNode));
			//		fT = Math.Atan((fY/fX) * Math.Sin(Trig.DegToRad(fPIncl)));
			fT = Math.Atan((fY/fX));
			fT = Trig.RadToDeg(fT) ;
			fT = Trig.TanQuadrant(fX, fY, fT);
			fL1 = fT + fPAscNode;
			fL1 = Trig.PutIn360Deg(fL1);
			fR1 = fPR * Math.Cos(Trig.DegToRad(fPsi));

			if (bInnerPlanet == true)
			{
				fA = Math.Atan((fR1 * Math.Sin(Trig.DegToRad(fEL - fL1)))/(fER - (fR1 * Math.Cos(Trig.DegToRad(fEL - fL1)))));
				fA = Trig.RadToDeg(fA);
				fLambda = 180 + fEL + fA;
			}
			else
			{
				fLambda = Math.Atan((fER * Math.Sin(Trig.DegToRad(fL1 - fEL)))/(fR1 - (fER * Math.Cos(Trig.DegToRad(fL1 - fEL)))));
				fLambda = Trig.RadToDeg(fLambda) + fL1;
				fLambda = Trig.PutIn360Deg(fLambda);
			}

			fPhase = 0.5 * (1 + Math.Abs(Math.Cos(Trig.DegToRad(fLambda - fPL))));

		}

		public static void CalcPlanetMag(DateTime dEpoch, DateTime dDate, double fPOrbPeriod, double fPEcc, double fPEpochLong, double fPPeriLong, double fPSMA, double fPAscNode, double fPIncl, double fEOrbPeriod, double fEEcc, double fEEpochLong, double fEPeriLong, double fESMA, bool bInnerPlanet,  double fMag1AU, ref double fMag)
		{
			double fD, fPN, fPM, fPL, fPV, fPR;
			double fEN, fEM, fEL, fEV, fER;
			double fPsi, fX, fY, fA, fL1, fR1, fLambda, fT, fPhase, fDistance;

			fD = UraniaTime.GetDaysBetween(dDate, dEpoch);
			fPN = (360.0/365.242191) * (fD / fPOrbPeriod);
			fPN = Trig.PutIn360Deg(fPN);
			fPM = fPN + fPEpochLong - fPPeriLong;
			fPL = fPN + ((360.0 / Math.PI) * fPEcc * Math.Sin(Trig.DegToRad(fPM))) + fPEpochLong;
			fPL = Trig.PutIn360Deg(fPL);
			fPV = fPL - fPPeriLong;
			fPR = (fPSMA * (1 - (fPEcc * fPEcc))) / (1 + (fPEcc * Math.Cos(Trig.DegToRad(fPV))));

			fEN = (360.0/365.242191) * (fD / fEOrbPeriod);
			fEN = Trig.PutIn360Deg(fEN);
			fEM = fEN + fEEpochLong - fEPeriLong;
			fEL = fEN + ((360.0 / Math.PI) * fEEcc * Math.Sin(Trig.DegToRad(fEM))) + fEEpochLong;
			fEL = Trig.PutIn360Deg(fEL);
			fEV = fEL - fEPeriLong;
			fER = ((1 - (fEEcc * fEEcc))) / (1 + (fEEcc * Math.Cos(Trig.DegToRad(fEV))));

			fPsi = Math.Asin(Math.Sin(Trig.DegToRad(fPL - fPAscNode)) * Math.Sin(Trig.DegToRad(fPIncl)));
			fPsi = Trig.RadToDeg(fPsi);

			fY = Math.Sin(Trig.DegToRad(fPL - fPAscNode));
			fX = Math.Cos(Trig.DegToRad(fPL - fPAscNode));
			//		fT = Math.Atan((fY/fX) * Math.Sin(Trig.DegToRad(fPIncl)));
			fT = Math.Atan((fY/fX));
			fT = Trig.RadToDeg(fT) ;
			fT = Trig.TanQuadrant(fX, fY, fT);
			fL1 = fT + fPAscNode;
			fL1 = Trig.PutIn360Deg(fL1);
			fR1 = fPR * Math.Cos(Trig.DegToRad(fPsi));

			if (bInnerPlanet == true)
			{
				fA = Math.Atan((fR1 * Math.Sin(Trig.DegToRad(fEL - fL1)))/(fER - (fR1 * Math.Cos(Trig.DegToRad(fEL - fL1)))));
				fA = Trig.RadToDeg(fA);
				fLambda = 180 + fEL + fA;
			}
			else
			{
				fLambda = Math.Atan((fER * Math.Sin(Trig.DegToRad(fL1 - fEL)))/(fR1 - (fER * Math.Cos(Trig.DegToRad(fL1 - fEL)))));
				fLambda = Trig.RadToDeg(fLambda) + fL1;
				fLambda = Trig.PutIn360Deg(fLambda);
			}

			fPhase = 0.5 * (1 + Math.Abs(Math.Cos(Trig.DegToRad(fLambda - fPL))));
			fDistance = Math.Sqrt((fER * fER) + (fPR * fPR) - (2.0 * fER * fPR * Math.Cos(Trig.DegToRad(fPL - fEL))));
	
			fMag = (5 * Math.Log10((fPR * fDistance)/(Math.Sqrt(fPhase)))) + fMag1AU;

		}


		public static void CalcPlanetDistance(DateTime dEpoch, DateTime dDate, double fPOrbPeriod, double fPEcc, double fPEpochLong, double fPPeriLong, double fPSMA, double fPAscNode, double fPIncl, double fEOrbPeriod, double fEEcc, double fEEpochLong, double fEPeriLong, double fESMA,  ref double fDistance)
		{
			double fD, fPN, fPM, fPL, fPV, fPR;
			double fEN, fEM, fEL, fEV, fER, fRho;
			
			fD = UraniaTime.GetDaysBetween(dDate, dEpoch);
			fPN = (360.0/365.242191) * (fD / fPOrbPeriod);
			fPN = Trig.PutIn360Deg(fPN);
			fPM = fPN + fPEpochLong - fPPeriLong;
			fPL = fPN + ((360.0 / Math.PI) * fPEcc * Math.Sin(Trig.DegToRad(fPM))) + fPEpochLong;
			fPL = Trig.PutIn360Deg(fPL);
			fPV = fPL - fPPeriLong;
			fPR = (fPSMA * (1 - (fPEcc * fPEcc))) / (1 + (fPEcc * Math.Cos(Trig.DegToRad(fPV))));

			fEN = (360.0/365.242191) * (fD / fEOrbPeriod);
			fEN = Trig.PutIn360Deg(fEN);
			fEM = fEN + fEEpochLong - fEPeriLong;
			fEL = fEN + ((360.0 / Math.PI) * fEEcc * Math.Sin(Trig.DegToRad(fEM))) + fEEpochLong;
			fEL = Trig.PutIn360Deg(fEL);
			fEV = fEL - fEPeriLong;
			fER = ((1 - (fEEcc * fEEcc))) / (1 + (fEEcc * Math.Cos(Trig.DegToRad(fEV))));
			
			fRho = Math.Sqrt((fER * fER) + (fPR * fPR) - (2.0 * fER * fPR * Math.Cos(Trig.DegToRad(fPL - fEL))));
			fDistance = fRho;
		}

		public static double CalcPlanetDiam(double fDistance, double fAngDiam1AU)
		{
			return (fAngDiam1AU / fDistance);

		}

		public static double CalcPlanetBrightLimb(double fSunRA, double fSunDecl, double fPRA, double fPDecl)
		{
			double fX, fY, fT, fDeltaRA;

			fSunRA = fSunRA * 15;
			fPRA = fPRA * 15;
			fDeltaRA = fSunRA - fPRA;

			fY = Math.Cos(Trig.DegToRad(fSunDecl)) * Math.Sin(Trig.DegToRad(fDeltaRA));
			fX = (Math.Cos(Trig.DegToRad(fPDecl)) * Math.Sin(Trig.DegToRad(fSunDecl))) - (Math.Sin(Trig.DegToRad(fPDecl)) * Math.Cos(Trig.DegToRad(fSunDecl)) * Math.Cos(Trig.DegToRad(fDeltaRA)));
			fT = Math.Atan(fY / fX);
			fT = Trig.RadToDeg(fT);
			return Trig.TanQuadrant(fX, fY, fT);
		}
	
	}
}
